#Tai-Hsien Ou Yang

library('mecdf')


#/nfs/apps/R/2.9.0/bin/Rscript

setwd("/home/to2232/likelihood")

FILE_INPUT_PED="gparchild.lst"
FILE_INPUT_NAME ="PopIndex.txt"
FILE_INPUT_LEN ="Numbers.txt"
FILE_OUTPUT_X ="x.gpn.txt"
FILE_OUTPUT_LK ="lk.gpn.txt"

parchild<-scan(FILE_INPUT_PED , list(id="",pr=""))
idx<-scan(FILE_INPUT_NAME , list(id=""))
ibdlen<-read.table(FILE_INPUT_LEN)

len12<-c()
len13<-c()

for(i in 1:length(idx$id))
{
templen=0

	pos1<-grep(parchild$id[i],idx$id)
	pos2<-grep(parchild$pr[i],idx$id)

	for(j in 1:length(pos2))
	{
		if(length(ibdlen[pos1[1],pos2[j]])>0&&is.na(ibdlen[pos1[1],pos2[j]])==0)
		{
			if(ibdlen[pos1[1],pos2[j]]>templen)
			{
				templen=ibdlen[pos1[1],pos2[j]]
			}
		
		}
	}

	len12=rbind(len12,templen)

templen=0
	pos1<-grep(parchild$id[i+1],idx$id)
	pos3<-grep(parchild$pr[i+1],idx$id)
	
	
	for(j in 1:length(pos3))
	{
		if(length(ibdlen[pos1[1],pos3[j]])>0&&is.na(ibdlen[pos1[1],pos3[j]])==0)
		{
			if(ibdlen[pos1[1],pos3[j]]>templen)
			{
				templen=ibdlen[pos1[1],pos3[j]]
			}
		}
	}
	
	len13=rbind(len13,templen)

}


x<-matrix(c(len12,len13), nrow = length(len12), ncol=2)

write.table(x, file = FILE_OUTPUT_X,col.names=F,row.names=F)

pcdf<-mecdf (x, continuous=FALSE, validate=TRUE,  project=FALSE, expandf=0.1)




lkmatrix<-list(id1=0,id2=0,id3=0,length=0)
for(i in 1:length(idx$id)) #
{
	for(j in 1:length(idx$id))
	{
		for(k in 1:length(idx$id))
		{
		
		lk<-pcdf(matrix( c(ibdlen[i,j],ibdlen[i,k]), ncol=2))
		lkmatrix<-rbind(lkmatrix,list(id1=i,id2=j,id3=k,length=lk))
		}
	}
}

write.table(lkmatrix, file = FILE_OUTPUT_LK,col.names=F,row.names=F)

